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Introduction 


In  the  following  we  consider  the  parameter-dependent  nonlinear  problem 


(1.1)  G(u,a)  =  0, 

where  G  :  Xm  x  +  Xm,  X  a  suitable  function  space  and  a  a  vector  of  real 
parameters.  (1.1)  will  typically  represent  a  parameter-dependent  nonlinear 
elliptic  system  of  dimension  m  and  order  two.  Important  examples  are  the  VLSI 
device  simulation  equations 


gi(u,v,w)  =  -  A  u  +  eu  v  -ew_u  -kj  =  0  , 

(1.2)  g2(u,v,w)  =  V*uneu-v  V  v  +  k2  =  0  , 

u 

g3(u,v,w)  =  -  V'u^e  V  w  +  k2  =  0  . 

In  this  case  the  parameters  a  enter  the  boundary  conditions  if,  for  example, 
current-voltage  characteristics  are  to  be  determined  (cf.[3]). 

The  solution  manifold  G  ^(0)  of  (1.1)  may  have  a  very  complex  structure. 

In  practical  applications  it  is  frequently  desirable  and  sufficient  to  compute 

certain  cross-sections  of  it,  i.e.  all  parameters  ot  are  kept  fixed  at  certain 

values  and  for  simplicity  suppressed  except  one  which  we  will  call  X.  A  simple 

graphical  representation  of  this  cross-section  may  be  obtained  if  a  functional 

% 

of  the  solution  is  depicted  versus  X.  It  should  be  noted,  however,  that  usually 
in  these  applications  such  functionals  are  in  fact  important  pieces  of 
information,  such  as  the  currents  for  (1.2). 


I 


We  assume  that  a  vector 


(1.3)  R(u(X)),  R  :  Xm  -*■  Rk 

of  differentiable  functionals  is  given  and  that  most  of  the  important 
information  about  the  problem  can  be  obtained  from  the  k  *  p  diagrams 

(1.4)  Ra  =  RjMc^)),  i  =  1,  ...»  k,  j  =  1,  ...,  p 
and  the  corresponding  solutions. 

The  solutions  u(X),  X  =  ot^,  in  (1.4)  and  from  them  the  curves  R^(u(X)) 
usually  have  to  be  computed  in  a  continuation  fashion.  Starting  from  a  known 
solution  u(X0)  another  solution  is  computed  by  utilizing  a  predictor-corrector 
scheme.  Frequently  an  initial  solution  to  start  this  continuation  procedure  is 
not  known  and  has  itself  to  be  determined  from  an  approximate  solution.  Along 
the  solution  curves  singular  points  have  to  be  expected  and  our  goal  is  to 
provide  methods  to  locate  these  points,  to  overcome  turning  points  and  to  switch 
branches  at  bifurcation  points. 

In  the  following  section  we  will  outline  a  continuation  algorithm  that  has 
been  successfully  used  for  the  scalar  case  (m  =  1  in  (1.1)).  It  provides  an 
efficient  and  robust  method  and  has  been  implemented  for  a  general  class  of 
second  order  elliptic  boundary  value  problems  discretized  by  finite  elements  in 
[2].  The  continuation  method  is  applied  on  a  coarse  triangulation  of  the  given 
domain  and  fine  grid  approximations  at  arbitrarily  specified  points  are  obtained 


2.  The  Continuation  Method 


In  the  following  we  assume  that  a  choice  of  the  parameter  X  =  a  in  (1.1), 

J 

(1.3)  and  also  of  the  functional  R^  in  (1.4)  has  been  made  and  denote  (1.1)  by 


(2.1)  G(u, X)  =0,  G  :  X™  x  R  +  Xm, 

and  by  r. 

The  parametrization  of  the  solution  u  of  (2.1)  by  X  and  thus  of  r  by  X  need 
not  be  possible  in  general.  Let  s  denote  the  arclength  along  a  solution  curve 
(u,X)  =  (u(s),  X(s))  of  (2.1).  Then  under  sufficient  smoothness  and  regularity 
assumptions  this  latter  dependence  is  differentiable  and  thus  (G°  =  Gu(uQ ,  X0 ) 
etc. ) . 


(2.2a)  Gj  u0  +  G®  X0  =  0, 

(2.2b)  !tu0"2  +  Xg  =  1 

in  a  solution  (uo,Xg)  =  (u(sg),  X(so))  of  (2.1).  Here,  subscripts  stand  for 
partial  derivatives  and  a  dot  for  differentiation  with  respect  to  s. 

In  singular  points  (uo,Xg)  of  (2.1)  the  linear  operator  G°  in  (2.2a)  is 
singular  and,  for  example,  a  solution  by  Newton's  method  breaks  down.  This 
singularity  may  be  removed  or  reduced  by  augmenting  (2.1)  by  a  normalizing 
equation,  a  well-known  example  of  which  is 


which  approximates  (2.2b).  o  is  called  the  pseudo-arclength  parameter.  6 
satisfying  0  £  0  £  2  may  be  used  as  a  weighting  parameter. 

Instead  of  (2.3)  we  propose  to  use 

(2.4)  Nr(u,X,o)  =  erQ(r-r0)  +  (2-0 >1Q ( X-XQ )  -  So  =  0, 

where  rQ  =  r(uQ),  rQ  =  r(u))a  _  and  r  denotes  the  above  functional.  From 

(2.4)  we  see  that  the  three  values  0  =  0,1,2  are  special  in  the  sense  that  for 

8  =  0(2)  the  augmenting  equation  characterizes  points  with  a  fixed  X(r)  -  value, 
while  for  0=1  this  point  lies  on  a  hyperplane  orthogonal  to  the  solution  arc 
in  (ro,X0). 

This  suggests  to  use  target  values  in  the  two  (physically  relevant) 
variables  r  and  X.  Such  a  strategy  seems  appropriate  as  long  as  starting  from  a 
point  on  the  solution  curve  another  point  corresponding  to  a  suitably  chosen 
target  value  can  be  computed  in  a  small  number  of  (corrector)  iterations.  This 
was  achieved  for  the  scalar  case  of  (2.1)  and  the  choice 

(2.5)  r(u)  =  lul. 

The  following  class  of  nonlinear  boundary  value  problems  in  the  plane 

-  V  •  a(x,y,u,Vu,X)  +  f(x,y,u,Vu,X)  =0  in  Q, 

(2.6)  u  =  gi(x,y,X)  on  W*  , 

a  •  n  =  g2(x,y,u,X)  on  3^2  =  3QS3Qi  , 

2 

was  considered.  Here,  Q  is  a  connected  domain  in  R  with  boundary  sq,  n  is  the 


V  J  ^.-4- # 


unit  normal  vector  on  3£1,  a  =  (ai,a2)  and  ai ,a2,f ,gi  ,g2  are  given  scalar 
functions. 

In  [2]  the  program  PLTMG,  an  implementation  of  an  adaptive  finite  element 
multi-grid  method  combined  with  the  continuation  technique  presented  below  is 
described  in  detail.  In  addition  to  some  theoretical  results  on  the 
continuation  method  [12,  13]  contain  numerical  results  for  several  problems  from 
the  applications  including  problems  with  turning-points,  symmetry-breaking  and 
other  bifurcation  points  as  well  as  parameter-switching  for  cases  depending  on 
several  (p  >  1  in  (1.1))  parameters.  [13]  also  contains  a  comparison  with  the 
classical  pseudo-arclength  continuation  method.  Other  problems  to  which  the 
program  has  been  applied  successfully  are  Steklov  eigenvalue  problems  for  which 
the  Dirichlet  condition  in  (2.6)  is  parameter-dependent. 

The  predictor  step  which  is  crucial  for  the  performance  of  the  method 
determines  predicted  values  of  the  form 


(2.7) 


(1+y)uo  +  Bu0 


Xq  +  aX0 


such  that  a,B, y  satisfy  the  three  equations 


N(u  ,X  ,a)  =  0, 
P  P 


(2.8) 


(u0,G(up,Xp) )  =  0, 


(ug ,G(up, Xp) )  =  0 


in  a  least-squares  sense.  Here  it  is  assumed  that  X  is  a  Hilbert  space  and 


)  denotes  its  inner  product. 

In  the  corrector  iteration  the  augmented  system 


(2.9) 


F(y)  = 


G(u,X) 

N(u,X,a) 


=  0,  yT  =  (uT,X) 


is  solved  by  Newton’s  method  starting  in  (u^X^).  The  restriction  of  the 
problem  (2.1)  to  the  (r,X)-plane  in  case  N  =  Nr  is  used  in  (2.8)  introduces  an 
additional  singularity  in  exceptional  points. 

If  in  a  solution  (uo,  X0)  simultaneously  tq  =  0  and  X0  =  0  then  the 
normalizing  condition  (2.4)  is  not  well-defined.  This  situation  occurs  in 
symmetry-breaking  pitchfork  bifurcation  points,  see  [11].  We  propose  to  switch 
temporarily  to  the  augmenting  equation  (2.3)  which,  as  explained  in  [12]  may  be 
viewed  as  asymptotically  equivalent  to  a  version  of  the  generalized  inverse 
iteration  ([10])  exploiting  symmetries  in  such  bifurcation  points.  The  precise 
choice  of  N  is  given  by  (cf.  [13]). 


(2.10) 


jro|  <  e  KugH  :  N  =  N^,  0*0, 
j r0 1  >  (1 -e )  "u0»  :  N  =  Nr,  9  *  0,  8  =  0, 


e  tuo"  £  r0  <  (1-e)  «u0"  :N=N,0£9<2. 
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3.  The  Corrector  Iteration 


A  point  (u^,  ^p)  computed  by  the  predictor  of  the  previous  section  is  used 
as  a  starting  guess  for  an  iterative  solution  of  the  augmented  system  (2.9). 
This  is  done  by  Newton's  method  and  thus  the  Dacobian 


has  to  be  evaluated  and  inverted. 

While  (2.10)  prescribes  which  N  to  choose  it  does  not  completely  define  the 
augmenting  equation.  The  weighting  parameter  A  still  has  to  be  picked  and  we 
assume  that  r  respectively  X  is  not  used  as  parameter  to  continue  past  a 
turning  point  in  r  respectively  X.  This  yields  9  *  0  if  X0  =  0  and  9  *  2  if 
ro  =  0.  Under  these  assumptions  it  was  shown  in  Proposition  4.2  of  [13]  that  F^ 
is  regular  except  in  bifurcation  points. 

When  computing  in  a  neighborhood  of  a  bifurcation  point  the  growth  of 
"F^-1lI  is  crucial  for  the  convergence  of  Newton  or  other  iterative  methods  used 
as  corrector,  cf.  [7],  We  give  here  the  detailed  proof  of  Proposition  4.3  in 
[13].  Let  locally  near  the  simple  bifurcation  point  yg  =  (uo,X0)  solutions  be 
parametrized  by  s  such  that  without  restriction  of  generality  yo  =  (u(0),X(0)). 
We  assume  in  the  following  that  yo  =  (u(0),X(0))  is  a  nondegenerate  ray  in  the 
sense  of  [8]. 

Theorem  3.1  If  the  normalizing  condition  N  is  chosen  according  to  (2.10)  then 
on  any  of  the  branches  passing  through  the  simple  bifurcation  point 
yo  =  (u(0),X(0))  of  (2.1)  for  G  >  0  sufficiently  small 


t.  ■»*#.  «1l m*rn.  ail 


""TW.KT.V*.  ^ 


(3.2) 


«Fy-1(y(s))l  =  0(  sj"1),  0  <  s|  <  6. 


Proof.  As  mentioned  in  [13]  we  have  to  generalize  the  proof  of  Theorem  1  in 
[8].  This  proof  covers  the  case  N  =  N  so  that  it  remains  to  consider 

<j 

N  =  N  .  We  note  that  the  row-vector  ao  =  (N  ,N.)  in  (3.1)  is  not  constant  in 
r  u  a 

this  case,  in  fact 


(3.3) 


a»(s>  ■  ■  <2-e»°> 


Let  us  further  denote  H(y)  =  G(u,X)  and  thus 


H'  =  (G  ,G.  )  :  Rn+1  +  Rn,  F  =  (H') 
u  ’  X  ’  y  ^  a0 ; 


As  in  [8]  we  introduce  the  one-dimensional  space  Xj  C  Rn+  spanned  by 
yo  =  (uo,Xq)^  on  one  of  the  branches,  the  one-dimensional  subspace  X2  c  Rn+^ 
such  that  the  kernel  K(H'(yo))  satisfies  K  =  X^  @  X2,  and  the  (n-1 )-dimensional 
subspace  X3  such  that  Rn+^  =  Xj  ©  X2  @  X3.  Similarly  let  Yj  =  R, 

Y3  =  R(H'(yo)),  where  R(A)  denotes  the  range  of  the  linear  operator  A,  and  let 
Y2  be  the  one-dimensional  subspace  of  Rn  such  that  Rn  =  Y2  ©  Y3 .  We  define  the 
orthogonal  projectors  Q  :  Rn  -*■  Y3  and  I  -  Q  :  Rn  +  Y2.  For  0  <  |s|  <5  we 
define  the  linear  operator  B(s)  :  Xi  x  X2  x  X3  Yj  x  Y2  x  Y3  by 

B(s) (xi  ,X2,X3)  =  (ao(s)(xi+X2+X3),  (I-Q )H ' (y( s) )  (xi +X2+X3  ) ,  QH ' (y(s) ) (xi +X2 +X3  ) ) . 


B(s)  is  equivalent  to  F^(y(s))  in  the  sense  that 
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where  C  and  D  are  regular  matrices  independent  of  s.  B(s)  can  be  decomposed  in 


a  natural  way  into  nine  suboperators  B^(s)  :  ♦  Y^,  \  <  i,j  <  3.  From  the 

above  definitions  it  follows  that  B2i(0),  B22(0),  831(0),  632(0)  and  823(0) 
vanish  identically.  From  (2.10)  we  conclude  the  crucial  fact 

•  *2  •  2 
ao(0)yo  =  9ro  +  (2-© )X0  *  0 


This  and  the  definition  of  X3,  Y3  imply  that  Bn(0)  and  833(0)  define 
bijections,  while  as  in  [8]  it  follows  that  82 y  (0 )  =  0  and  that  ^(O)  is  a 
bijection  between  X2  and  Y2.  We  can  thus  expand  B(s)  in  a  neighborhood  of  s  =  0 
as 


B(s) 


Bu(0) 

2B2i(0) 

sB31(0) 


B12(0) 

sB22(0) 

SB32 (0) 


+  h.o.t. 


where  B..(0),  B..(0)  and  8. .(0)  are  linear  operators  from  X.  into  Y.  which  are 
xj  ’  ij  K  i  j 

uniformly  bounded  with  respect  to  S.  Thus  det  B(s)  =  0(|s|)  and  hence  B'1 (s) 
can  be  written  as 


/En(s) 

B-1(s)  =  f  sE2 1 (s) 
\  SE3 1 (s) 


s_1E12(s)  Ei 3 (s)\ 

s-1E22(5)  E2 3 ( s )  j 

E  3  2 (s )  E33 (s)/ 


where  similarly  the  E^j  are  uniformly  bounded  and  E22(0)  =  (822(0) )-1  is 

regular. 


4.  Spurious  Solutions 

It  Is  well-known  that  discretized  problems  possess  frequently  many  more 
solutions  than  the  continuous  problems.  All  discrete  solutions  that  do  not 
correspond  to  solutions  of  the  continuous  problem  and  disappear  for  h  ♦  0,  h  the 
discretization  parameter,  are  called  spurious  or  numerically  irrelevant 
solutions.  If  we  consider  the  simple  model  equation 

(4.1)  -  Au  =  Xf(u) 

then  there  are  essentially  two  types  of  those  solutions.  One  type  bifurcating 
from  X  =  ±  <*>  in  case  f(u)  has  zeros  and  the  other  bifurcating  from  Hun  =  “  if 
the  growth  of  f  is  sufficiently  strong.  It  might  be  argued  that  not  too  much 
effort  should  be  spent  on  investigating  those  spurious  solutions  if  the  main 
interest  is  that  of  approximately  solving  the  continuous  problem.  There  are, 
however,  at  least  two  counter  arguments. 

First,  in  a  multi-grid  context  where  one  hopes  to  be  able  to  use  rather 
crude  coarse  grids  on  which  the  continuation  will  be  done,  one  has  to  be  aware 
of  the  spurious  solutions  since  for  'large'  h  they  are  present  for  'small' 
values  of  *ufl  respectively  X  and  thus  may  be  more  easily  encountered.  Second, 
it  has  been  shown  recently  ([9])  that  for  the  second  type  of  spurious  solutions 
some  of  these  may  bifurcate  from  a  relevant  branch.  Thus  they  definitely  will 
be  encountered  and  it  becomes  necessary  to  distinguish  between  spurious  and 
relevant  solutions. 


If  we  specialize  (4.1)  to 


then  the  proof  in  [9]  covers  the  case  n  =  1  and  a  uniform  grid  of  an  even  number 
of  points.  The  midpoint  of  [0,1]  where  the  relevant  solution  attains  its 
maximum  is  then  not  a  grid  point.  The  discrete  solution  attains  its  maximum  in 
the  two  neighboring  points  and  from  this  branch  bifurcate  two  branches  attaining 
their  maximum  in  either  of  these  points. 

We  had  observed  these  solutions  and  their  bifurcation  previously  when 
numerically  solving  (4.2)  for  n  =  2  by  discretizing  with  the  standard  difference 
method  on  a  corresponding  grid.  We  also  have  applied  a  multi-grid  method  to 
compute  spurious  solutions  ([14]).  Since  our  goal  is  not  to  contribute  to  the 
knowledge  about  the  spurious  solutions  but  rather  to  provide  a  general  purpose 
and  if  possible  foolproof  program  to  solve  a  general  class  of  problems,  it  is 
natural  to  investigate  the  above  situation  with  the  program  PLTMG. 

Our  expectation  is  that  the  implemented  continuation  method  allows  us  to 
follow  the  spurious  solution  branches  because  these  are  genuine  solutions  of  the 
discrete  problems.  On  the  other  hand  the  multi-grid  process  using  several 
discretizations  simultaneously  could  help  to  identify  spurious  solutions.  But 
additionally  all  the  other  features  offered  by  the  program  as  the  adaptive 
refinement  and  the  error  estimators  will  have  their  impact  on  the  solution 
process  respectively  may  give  other  hints  about  which  solutions  are  relevant. 

Since  the  program  uses  a  finite  element  discretization  with  piecewise 
linear  elements  we  try  to  come  close  to  the  above  discretization.  The  unit 
square  is  subdivided  into  nine  subsquares  of  length  h  =  1/3  and  then  all  these 
squares  are  divided  into  two  triangles  by  the  SW  -  NE  diagonal.  This  way  the 
Laplace  operator  will  be  discretized  in  the  same  way  as  by  the  five  point  star 
difference  method.  The  right-hand  side,  however,  in  the  finite  element  context 
discretized  by  a  quadrature  formula  will  be  different  from  the  pointwise 


discretization  of  the  difference  method.  This  will  introduce  a  bias  in  the 
SW-NE  direction  and  we  do  not  expect  the  relevant  solution  to  attain  its  maximum 
in  the  four  interior  nodes. 

While  the  numerical  results  presented  in  the  following,  of  course,  make  no 
statements  about  the  mathematical  properties  of  the  problem  they  provide  some 
valuable  insight  into  the  practical  use  of  the  program  and  the  numerical 
solution  in  the  presence  of  spurious  solutions. 

The  branch  starting  at  A  =  0,  r  =  Hu*  =  0  for  the  above  discretization  with 
four  nodes  attains  its  maximum,  as  expected,  in  the  points  (1/3,  1/3),  (2/3, 

2/3)  and  has  two  equal  but  smaller  values  at  (1/3,  2/3),  (2/3,  1/3).  When 
following  this  curve  beyond  the  limit  point  another  sign-change  of  det  (G  )  is 
encountered.  Even  a  tracing  of  the  neighborhood  of  this  point  with  small 
continuation  steps  shows  an  abrupt  change  of  the  sign.  When  asked  to  locate 
this  point  the  program  ends  up  on  a  different  branch  and  even  a  subsequent 
branch  switching  does  not  lead  back  to  the  relevant  branch.  This  suggests  that 
instead  of  bifurcation  a  perturbed  bifurcation  seems  to  be  present.  This  is  not 
so  relevant  here,  since  the  fact  remains  that  a  sign-change  seems  to  indicate 
bifurcation  and  naturally  a  spurious  branch  is  found.  Table  4.1  shows  the 
continuation  to  and  the  location  of  the  'bifurcation'  point.  In  the  tables  NI 
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denotes  the  number  of  corrector  iterations  and  is  a  quantity  that  changes  sign 
at  simple  bifurcation  points  (cf.  [2]). 


In  the  following  we  present  multi-grid  results  for  four  different  cases. 
First  we  have  followed  the  relevant  branch  up  to  r  =  4.5  respectively  r  =  5  and 
computed  a  3-grid  solution  for  fixed  norm  with  uniform  and  adaptive  refinement. 
Table  4.2  contains  only  the  refinement  results. 


N1 

X 

Hu" 

X 

• 

r 

det(G  ) 
u 

6 

Continue  to  X  = 

.1  then  to  r 

=  1,  r  =  3, 

r  =  3.44,  r 

=  3.46. 

1 

0.000E+00 

0.000E+00 

0.995E+00 

0.000E+00 

0.492E 

5 

0.200E+01 

1 

0.100E+01 

0.11 6E+00 

0.992E+00 

0.128E+00 

0.276E 

5 

0.180E+01 

1 

0.390E+01 

0. 100E+01 

0.244E+00 

0.970E+00 

0.190E 

3 

0. 129E+00 

2 

0.156E+01 

0.300E+01 

-0.756E+00 

0.653E+00 

-0.62 3E 

3 

0.592E+00 

2 

0.111E+01 

0.344E+01 

-0.669E+00 

0.741E+00 

-0.1 85E 

2 

0.1 36E-01 

1 

0.109E+01 

0.346E+01 

-0.664E+00 

0.745E+00 

0.176E 

2 

-0.130E-01 

Determinant  changes  sign;  find  bifurcation  point  by  secant  method  on  6. 

1 

0.110E+01 

0.345E+01 

-0.121E-01 

0.130E-01 

-0.113E 

-1 

0.809E-05 

1 

0. 1 10E+01 

0. 345E+01 

0. 192E+00 

-0.215E+00 

0.123E 

0 

-0.883E-04 

1 

0.1  IDE +01 

0.345E+01 

0.497E-06 

-0.591E-03 

0.482E 

-2 

-0.369E-09 

Table  4.1  Continuation  along  relevant  branch,  detection  of  spurious  branch. 

A  posteriori  error  estimates  based  on  solving  local  problems  with  higher-order 
elements  [5]  are  evaluated  by  the  program  and  in  Table  4.4  we  give  the  number  of 
estimated  correct  digits  in  three  different  norms. 

Following  the  irrelevant  branch  that  attains  its  maximum  at  (1/3,  1/3)  to 
r  =  5.5  and  again  refining  uniformly  and  adaptively  yielded  the  results  in  Table 
4.3.  We  see  that  for  the  adaptive  refinement  it  seems  impossible  to  conclude 
from  the  error  estimations  alone  that  the  solution  is  irrelevant,  while  for  the 
uniform  refinement  the  error  was  estimated  as  bigger  than  the  solution  resulting 
in  a  negative  number  of  estimated  correct  digits.  In  this  latter  case  also  the 
solution  on  the  finest  grid  did  not  attain  its  maximum  at  (1/3,  1/3)  but  closer 
to  the  center  of  the  domain;  i.e.  the  refinement  had  led  onto  a  different 
spurious  branch. 
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Level  1  solution  at  r  =  5. 

2 

0.264E+00 

0, 500E+01 

-0.255E+00  0. 962E+00 

0.865E 

4 

0.891E+00 

Level 

Refine  uniformly. 

2  solution,  no.  of  vertices 

49. 

6 

0.296E-01 

0. 500E+01 

-0. 566E-01  0.993E+00 

0.530E 

5 

0.885E+00 

Level 

3  solution,  no.  of  vertices 

169. 

8 

0.836E-02 

0.500E+01 

-0.230E-01  0.990E+00 

0.626E 

5 

0.251E+00 

Level  1  solution  at  r  =  4.5. 

2 

0.431E+00 

0.450E+01 

-0. 379E+00  0.921E+00 

0.393E 

4 

0.824E+00 

Level 

Refine  adaptively. 

2  solution,  no.  of  vertices 

51. 

5 

0.598E-01 

0.450E+01 

-0. 1 15E+00  0.983E+00 

0.334E 

5 

0.465E+00 

Level 

3  solution,  no.  of  vertices 

265. 

9 

0.402E-02 

0.450E+01 

-0.315E-01  0.862E+00 

0.482E 

5 

0. 504E-03 

Table  4.2  Results  of  refinement  on  relevant  branch. 

One  of  the  possible  ways  provided  by  the  program  to  depict  the  results 
graphically  is  to  generate  a  histogram.  Over  a  square  grid  numbers  from  0  to  9 
are  printed  covering  uniformly  the  range  of  the  solution  values  from  minimum  to 
maximum.  In  Figures4.1  and  4.2  we  see  which  solutions  were  computed  on  the 
coarsest  and  on  the  finest  grid  for  uniform  and  adaptive  refinement. 


NI  X  lull  X  r  det(Gu)  fi 


Level  1  solution  at  r  =  5.5. 


2 

0.560E-01 

0.550E+01 

-0.772E-01  0.994E+00 

-0.839E 

5 

0.219E+01 

Level  2 

Refine  uniformly. 

solution,  no.  of  vertices 

49. 

9 

0.328E-01 

0.550E+01 

-0. 367E-01  0. 508E+00 

-0.361E 

5 

-0. 1 17E+00 

Level  3 

solution,  no.  of  vertices 

169. 

10 

0.474E-02 

0.550E+01 

0. 1 12E-01  -0.990E+00 

0.615E 

5 

0.182E+00 

Level  2 

Refine  adaptively. 

solution,  no.  of  vertices 

41. 

4 

0.968E-02 

0.550E+01 

-0. 194E-01  0.994E+00 

-0.824E 

5 

0.691E+00 

Level  3 

solution,  no.  of  verticies 

174. 

7 

0.111E-03 

0.550E+01 

-0. 386E-03  0.999E+00 

0.638E 

5 

0.776E-01 

Table  4.3  Results  of  refinement  on  spurious  branch. 


Table  4.4  ■  Numbers  of  estimated  correct  digits. 


Norm 

Relevant 

Solution 

Spurious  Solution 

Uniform 

Adaptive 

Uniform 

Adaptive 

H1 

0.526 

0.593 

■■ 

0.448 

L2 

1.381 

1.268 

1.198 

L°° 

0.965 

1.186 

-1.099 

- - 

1.183 
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Figure  4.1  Histograms  of  the  spurious  solution  on  the  coarsest  grid 
and  for  uniform  and  adaptive  refinement. 
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Figure  4.2  Histograms  of  the  relevant  solution  on  the  coarsest  grid  and 
for  uniform  and  adaptive  refinement. 
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